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ABSTRACT 


Progress in research into fundamental mechanisms of combustion in high 
speed air flows is described. Three major areas are covered: the application 
of a modular combustor model for the liquid fueled sudden expansion combustor 
configuration to the analysis of flame stability and the interpretation of 
combustor test data; the extension of the modular model to the analysis of 
boron slurry fueled ramjet combustors; and the development of analytical models 
of the ducted rocket combustor. The modular approach is demonstrated to be 
a useful tool for the analysis and interpretation of sudden-expansion combustor 
flame stabilization phenomena and to be of considerable utility in the analysis 
of combustor test data. Models for the combustion of boron particles are 
reviewed and put into a form suitable for incorporation in a modular model of 
a boron slurry fueled combustor. Particle tracking methods suitable for use 
in the context of a modular model formulation have been defined and developed 
for the required application. The approach required to provide an analytical 
model of the three-dimensional recirculating flow in a ducted rocket combustor 
has been defined and the development of an aerodynamic model has been initiated. 






1. INTRODUCTION 


Design trends for current and future airbreathing propulsion systems are 
trending toward lower weight and smaller volume designs. These requirements, 
and increased demands for higher performance from airbreathing propulsion 
systems in general and ramjet devices in particular, have resulted in the need 
to upgrade existing technologies to meet the emergence of more stringent 
design requirements. Achieving reduced volume and weight requires the use of 
shorter combustors and high energy, high density fuels. Systems constraints 
have dictated the use of sudden-expansion (dump) combustors capable of 
operating effectively at high combustion intensities over wide ranges of 
conditions. Problems of flame stabilization, flame propagation and spray 
combustion have been encountered. To aid in the development of solutions to 
these problems, analytical combustor models are a requirement: the ability 
to compute in some detail combustion chamber flowfields is necessary in order 
to understand the phenomena that occur in existing combustors and to predict 
the performance of new combustor concepts. The insight gained through the use 
of analytical combustor models can be of substantial value in the planning of 
a combustor test program and in the interpretation of combustor and combustor 
component test data. 

In response to the needs just outlined, a detailed model of the sudden- 
expansion liquid fueled ramjet combustor has been developed, as part of earlier 
AFOSR-sponsored work (Ref. 1). This model makes use of the modular concept, 
in which the combustor flowfield, represented schematically in Fig. 1, is 
broken down into three major components: a directed flow, which is treated 
as parabolic, a recirculation zone, assumed to be represented by well- 
stirred reactor(s), and a turbulent shear layer along the dividing streamline 
which separates the other two regions. The shear layer serves as the coupling 
region between the other two model components; fluxes of species and energy 
across this shear layer form the boundary conditions on the two computational 
regions. Finite-rate chemistry, based on the quasiglobal model (Ref. 2) is 
included in the formulations for both the directed flow and well-stirred 
reactor regions, although for the modular model calculations described in this 
paper, the recirculation region well-stirred reactor formulation has been 
restricted to a global finite-rate chemistry model. The directed flow is 
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Figure 1. Schematic of Sudden-Expansion (Dump) Burner. 
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assumed to be fully turbulent, with the turbulent viscosity defined by a two- 
equation turbulence model (Ref. 3). A key feature of this approach is the 
provision for the shear layer coupling region in the model. Through the use 
of this element of the model, the division of the mass flux between the 
directed flow and the recirculation region is computed iteratively rather than 
specified empirically. Furthermore, the directed flow region is computed in 
detail as a two-dimensional parabolic flowfield, rather than through a one¬ 
dimensional approximation, allowing the use of detailed computations of the 
mixing and chemical reactions in this region of the combustor. 

A variety of results obtained in the process of validating the modular 
model approach for the prediction of ramjet combustor performance and flow- 
field characteristics have been described in Refs. 1 and 4. During the 
current program, attention has been focused on further development of certain 
of the fundamental elements of the model, with particular emphasis on the 
stirred-reactor representation of flame-stabilizing recirculation regions. 
Other aspects of the current work include the development of boron combustion 
models and the definition of modeling approaches for ducted rocket combustor 
configurations. 

Both gas generator ramjets (ducted rockets) and slurry-fueled ramjets 
are examples of propulsion devices which involve the turbulent mixing and 
combustion of multi-phase flows. In both of these applications boron-loaded 
fuels are of considerable interest because of their high energy content, and 
for both of these systems it is desirable to develop models which can aid 
in obtaining a fundamental understanding of the phenomena involved in two- 
phase mixing and combustion and which can be used to parametrically assess 
combustor performance. This need is being addressed in two of the tasks that 
are part of the current program. For boron slurry-fueled ramjets, the 
approach being taken involves a direct extension of the modular modeling 
approach developed for liquid-fueled ramjets and described in Refs. 1 and 4, 
while the application of modular modeling to the more complex geometry associ¬ 
ated with the ducted rocket configuration is also being investigated. In 
the latter work, the approach being followed involves the development of a 
three-dimensional elliptic Navier-Stokes solution technique to be used to 
delineate the flowfield regions that are to be described by separate elements 
of the modular approach. This solution procedure is being developed to 
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obtain the cold-flow aerodynamics of the combustor configuration; the detailed 
modeling of this configuration incorporating the effects of multiphase flow 
and chemical heat release is to be carried out using the modular approach. 

The three subsequent sections of this report reflect the three main 
areas of investigation outlined above. Thus, in Section 2, further examination 
of the use of the modular model of a liquid-fueled sudden expansion combustor 
is described, with emphasis on the prediction of flame stabilization phenomena. 
A detailed model of the combustion dynamics inherent in a boron slurry-fueled 
ramjet combustor requires the development of models for the combustion of boron 
particles and of techniques for treating the details of the combustion of 
boron particle clouds, and progress in this area is discussed in Section 3. 
Then, in Section 4, the technique being utilized to develop a gas-generator 
fueled ramjet (ducted rocket) aerodynamic model is described, and the overall 
progress made during this program is summarized in Section 5. In each of these 
areas close coordination has been maintained with government laboratory experi¬ 
mental work, including liquid-fueled ramjet and ducted rocket investigations 
underway at the Aeropropulsion Laboratory of the Air Force Wright Aeronautical 
Laboratory (AFWAL) and boron-slurry ramjet and ducted rocket work underway 
at the Naval Weapons Center (NWC) under joint AFOSR/NWC sponsorship. 




2. THE APPLICATION OF MODULAR MODELING TO RAMJET 


PERFORMANCE INVESTIGATIONS* 


The well-stirred reactor, shown schematics!ly in Fig. 2a, is a laboratory 
device in which very high mixing rates are achieved. In general, laboratory 
stirred reactors are designed to ensure that the mixture within the reactor 
is spatially uniform, so that unreacted feed material is continuously and 
uniformly mixed with combustion products, reacting for a time defined by the 
average residence time of the reactor before exiting. For the limit of perfect 
stirring, this average residence time is given simply by the ratio of the 
stirred reactor volume to the feed mass flux, V/m; the product mass flux is 
of course equal to the feed mass flux. In this limit the governing equations 
for the stirred reactor state reduce to algebraic relations (Ref. 4), which 
allows the use of rapid and efficient solution procedures, incorporating either 
global finite-rate kinetics or the more detailed quasiglobal kinetics formulation. 

In the lower half of Fig. 2, a typical recirculation region in a sudden- 
expansion combustor is sketched. Because of the high turbulence intensity and 
large mixing rates generated in the flow reversal within the recirculation 
the limiting behavior of this region can be thought of as well-mixed; that is, 
the state of the recirculation region is defined by the chemical kinetics of the 
reactions occurring and not by the mixing rate. In this sense the volume within 
the recirculation region sketched in Fig. 2b is similar to the volume within the 
laboratory stirred reactor of Fig. 2a, and the same solution technique can be 
used to obtain the thermochemical state. A comparison of Figs. 2a and 2b also 
shows the major difference between the recirculation region and the laboratory 
stirred reactor. In the laboratory stirred reactor discrete reactant and prod¬ 
uct streams can be identified, but in the well-stirred reactor model of a 
sudden-expansion recirculation zone, separate reactant and product streams cannot 
be defined. Instead, reactants enter the recirculation region and products leave 
the recirculation region by turbulent diffusion through the shear layer, inte¬ 
grated over the surface area of the dividing streamline, and the product stream 
is defined similarly. These fluxes are of course equal and opposite. 


* The material incorporated in this section is taken from "Interpretation of 
Ramjet Combustor Test Data," by P. T. Harsha and R. B. Edelman, AIAA Paper 
81-1433, to be presented at the AIAA/SAE/ASME 17th Joint Propulsion Conference, 
July 27-29, 1981. 
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(a) Schematic design of spherical 
well-stirred reactor. 



(b) Schematic of sudden-expansion 
recirculation region. 


FIGURE 2. Comparison of well-stirred reactor 
with recirculation region flowfield. 
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The fact that in the well-stirred reactor model of a sudden-expansion 
recirculation region the feed rates are defined by the fluxes of reactants 
and products through the shear layer region is the reason that this region 
has been identified as a key element of the modular model described in Refs. 

1-3 and 6. Nevertheless, this shear layer is modeled fairly crudely as a region 
of linear gradients: that is, the gradient in a quantity at the dividing 
streamline is approximated by 



where the subscripts P, R, and W refer to conditions in the directed flow 
immediately outside of the recirculation zone, in the recirculation zone itself, 
and along the dividing streamline, respectively. Moreover, the width of the 
shear layer is assumed to be given by the linear expression 

£ = ax + b (2) 

in which a and b are constants to be defined through comparison of model 
predictions with experimental data. 

As is described in more detail in Refs. 1, 4-6, the overall flowfield 
computation using the modular model proceeds as follows: a dividing stream¬ 
line shape is assumed, and the shear layer width expression and shear stress 
distribution along the dividing streamline are specified, the latter through 
use of a specified "skin friction" coefficient. An initial state for the 
stirred reactor is assumed, which defines an initial guess for the species 
distribution and temperature within the recirculation zone. With these 
available as boundary conditions, a finite-difference computation of the mixing 
region external to the recirculation is carried out to the end of the recircu¬ 
lation region. This calculation defines the species mass fraction and 
temperature gradients along the dividing streamline, and these values are used 
as outlined above to obtain the stirred reactor feed rates. A new stirred 
reactor computation is carried out using these feed rates, resulting in a new 
specification of the species mass fractions and temperature within the recircu¬ 
lation region, and the parabolic mixing computation is repeated. This procedure 
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is continued until changes in the stirred reactor composition and temperature 
from computation to computation become small, typically 0.1%, at which point 
the coupling iteration has converged. The parabolic calculation is then 
carried out to the end of the combustion chamber, completing the solution. 

The entrainment rate into the recirculation region is an important 
parameter in the overall combustor flowfield computation, and as the previous 
paragraph indicates, this parameter is not specified, but is calculated as a 
part of the modular model solution procedure. Thus it is of considerable 
interest to compare the entrainment rate computed by the model with experimental 
data. A variety of data for recirculation region entrainment rates was 
reviewed by Curran (Ref. 7) and the range of the observed correlations is 
shown in Fig. 3: these correlations are for both cold flow and hot flow. 

Also shown in Fig. 3 are the entrainment mass flux ratios predicted using the 
modular model in several different calculations, involving different area 
ratios and cold flow as well as reacting flow. Other results from many of 
these computations are described in Ref. 4. The range of correlations shown 
did not include data from the investigations for which the calculations shown 
in Fig. 3 were made; nevertheless, the predictions show reasonably good agree¬ 
ment with the observed correlations. Of particular significance is the fact 
that the predicted reduction of the entrainment rate with combustion compared 
to cold flow at the same area ratio is in agreement with the observations 
reported by Curran. This represents an important result, both from the 
standpoint of validation of the modular model approach, even though a rela¬ 
tively crude shear layer model has been used, and because it suggests a more 
systematic study that can be pursued using the analytical model to examine 
a wider range of operating conditions than those on which the correlations 
were originally based. 

One of the most important determinants of liquid fueled ramjet performance 
is the fuel distribution at the combustor entrance. For the sudden-expansion 
liquid fueled ramjet, the interaction between inlet fuel distribution and 
recirculation zone state is a further important consideration. A variety 
of observations have shown that for wall injection of the fuel, the recircula¬ 
tion zone equivalence ratio is generally more fuel-rich than the overall 
(global) equivalence ratio. This is demonstrated by the data of Schmotolocha 
and Economos (Ref. 8), which was obtained for three different inlet conditions 
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IGURE 3. Comparison of measured and 

nredicted entrainment rates. 







and at two different locations of the fuel injection orifices upstream of the 
dump plane. Indeed, the measurements described in Ref. 8 indicate that for 
overall equivalence ratios greater than 0.2, the equivalence ratio in the 
recirculation region is always greater than unity. These data are shown in 
Fig. 4, in which the results of a modular model computation at one overall 
equivalence ratio are also shown. As was the case for the entrainment rate 
pedictions, the agreement between the experimental results and the computation 
is reasonably good. 

In the experiments described in Ref. 8 plain orifice fuel injection was 
used. For plain orifices, fuel penetration characteristics are reasonably 
well documented, and an empirical penetration correlation was used to establish 
initial conditions for the computations used to obtain the result shown in 
Fig. 4. However, plain orifices are not widely used in ramjet design because 
of their poor turndown characteristics which result in the need for several 
sets of orifices to obtain a wide ramjet operating range. The most commonly 
used fuel injector in liquid fuel ramjet applications is the poppet injector, 
which has a much broader turndown range than does a plain orifice. However, 
penetration characteristics from poppet injectors under ramjet operating 
conditions are not well documented. Indeed, there exists one investigation 
of poppet (and other nozzle) penetration characteristics under ramjet operating 
conditions in which it was observed that the penetration from poppets is 
essentially nil. Ref. 9. While other results indicate that penetration 
characteristics similar to orifice injection can be obtained using poppets. 

Ref. 10, it remains clear that the penetration characteristics of poppets 
are not as well known as for plain orifices. 

At one operating condition, which involves a relatively high inlet 
temperature, poor ramjet combustor performance was observed in a recent 
development program: at this one condition performance was considerably 
poorer than at other test conditions. Review of the available results suggested 
that fuel distribution effects were a possible cause, and in particular that 
at the operating condition in question the fuel penetration from the wall- 
mounted fuel injectors may have been negligible. To investigate this possibility, 
and to study the effects of fuel distribution on combustor performance, a series 
of parametric modular model computations was carried out. 
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In the modular approach, fuel injection is modeled using an empirical 
correlation for penetration and breakup time correlation. The latter correla¬ 
tion defines the downstream position at which the penetration is computed; 
from this point the fuel spray is assumed to vaporize at a rate given by a 
bulk spray vaporization correlation. After injection, the fuel is assumed 
to reside in an annular region at the position specified by the penetration 
correlation, and spreading of this annulus is computed using a turbulent mixing 
model. For the conditions of interest in the parametric study, vaporization 
occurs rapidly, so that the parameter to be varied was the fuel penetration 
height. 

Fig. 5 shows the predicted recirculation zone equivalence ratio as a 
function of fuel penetration height obtained from the parametric study. Of 
particular interest is the result for near-zero penetration height, t/r ~ 0 . 
Under this condition, for a global equivalence ratio of unity, the predicted 
recirculation zone fuel/air ratio is extremely fuel-rich, with 4 ^ = 4.5. 

These high values of equivalence ratio produce a substantially cooler recircula¬ 
tion region than would be obtained for conditions nearer stoichiometric, and 
this in turn leads to reduced flame propagation rates and lower overall 
combustion efficiency in the combustor. Evaluation of the effect of the rich 
fuel/air ratio in the recirculation zone on overall performance is, however, 
complicated by difficulties that are encountered in modeling chemical kinetics 
in very fuel-rich regions. 

There is no currently available finite-rate chemical kinetics model which 
can provide adequate results for complex hydrocarbon fuels under fuel-rich 
conditions, at temperatures at which the effects of dissociation manifest 
themselves. The modular model formulation used for the computations described 
in this paper includes as options either a one-step finite-rate model (Ref. 4) 
or the quasiglobal kinetics model (Refs. 2, 14). In the former case, the 
products of combustion are the fully-reacted species Wfl and COg, so that the 
effects of dissociation are ignored, and the overall reaction rate is under¬ 
predicted. The quasiglobal model used in the modular formulation is valid 
for complex hydrocarbon fuels and includes the effects of dissociation, but 
is not appropriate for equivalence ratios much above stoichiometric. Other 
models have been proposed for fuel-rich conditions, such as the four-step 
model of Dryer and Glassman, Ref. 12, but this model ignores dissociation 
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FIGURE 5. Recirculation Zone Equivalence Ratio as a Function of 
Penetration Height 
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effects. Thus it is difficult to obtain reliable estimates for the chemical 
kinetic rates that pertain in combustors in which significant regions of the 
flowfield are highly fuel-rich. 

This difficulty is reflected in the performance results shown in Fig. 6. 
While the trend of combustion efficiency as a function of fuel penetration 
distance obtained through the use of the one-step global kinetics model is 
as expected, the magnitude of the overall change is relatively small and 
probably incorrect. The reason that the magnitude of these results can be 
expected to be incorrect is that the global finite-rate model can be expected 
to significantly overpredict the recirculation zone temperature for fuel-rich 
conditions, since the only effect of fuel-rich operation that this model 
recognizes is a dilution of the products of combustion with excess fuel for 
<t> » 1. With this model, the reaction rate is linear with respect to fuel 
concentration, whereas evidence exists that for fuel-rich conditions reaction 
rates do not continue to increase with fuel concentration, and, as noted 
earlier, dissociation effects are also ignored. 

Nevertheless these results indicate that for the conditions considered, 
low values of fuel penetration distance are directly coupled to reductions 
in overall combustor performance. Since further development testing using 
fuel injection schemes designed to increase fuel penetration resulted in 
substantial performance improvements, it can be considered that the trends 
demonstrated by this modeling effort have been experimentally verified. 

In the descriptions of the modular approach to the development of combustor 
models outlined in Refs. 1, 4, 6 and 7 it has been stressed that one of the 
features of the modular approach is that each of the elements involved in the 
overall combustor model can be independently developed. Just such an approach 
is being followed at SAI where, under a DOE-sponsored program, the extension 
of the quasiglobal model to fuel-rich conditions is being explored using the 
well-stirred reactor element of the model. One result of this work is shown 
in Fig. 7, on which predictions using several different kinetics models 
are compared to measured temperature data obtained in a laboratory well- 
stirred reactor burning toluene under both fuel-lean and fuel-rich conditions. 
All of these results are for the same stirred reactor residence time. 

For the fuel-rich condition. Fig. 7 shows that the extended quasiglobal model 
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FIGURE 6. Predicted overall combustion efficiency as a 
function of fuel penetration height. 








substantially improves the stirred reactor temperature prediction as compared 
to either the original quasiglobal formulation (Refs, 2, 11) or the four-step 
model (Ref. 12). 


Although the extended quasiglobal formulation has not yet been fully 
developed, an indication of the effect of using this approach in a combustor 
calculation is demonstrated by Fig. 8, which shows recirculation zone tempera¬ 
tures for the conditions of the preceding section, computed using the four-step 
model (Ref. 12), the one-step model (Ref. 4), and the extended quasiglobal 
formulation (Ref. 13). While these results are provisional in that the 
extended quasiglobal model has not at this point been verified for ramjet fuels 
or operating conditions, they indicate a significant effect on recirculation 
zone temperature. In particular, at 4>^ z = 3.0, the recirculation zone tempera¬ 
ture predicted using the extended quasiglobal model is some 300°K lower than 
that obtained from the one-step model. This difference, coupled with the 
radically different species distributions obtained from the two models, can 
be expected to exert a significant influence on predicted flame propagation 
rates. 


In summary, a comparison of modular model predictions with available 
data shows that the model is capable of predicting both the entrainment rate 
into sudden-expansion combustor recirculation regions and the overall fuel/air 
ratio in these regions. Both of these parameters are crucial to investigations 
of flame stability in a liquid fueled sudden-expansion combustor. An appli¬ 
cation of the approach in the interpretation of ramjet combustor test data 
has also been demonstrated: This demonstration shows that the model provides 
a means of verifying suspected causes of combustor performance deficiencies. 

The inadequacy of current chemical kinetics models under very fuel-rich condi¬ 
tions in circumstances under which the effects of dissociation are apparent 
provides a current limitation on the use of the modular model approach. 

However, related chemical kinetic model development work which can remove this 
limitation has been described. 
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3. MODIFICATIONS TO MODULAR MODEL FOR BORON SLURRY FUEL COMBUSTION 


Boron slurry fuels have a high energy density and are therefore extremely 
attractive for volume-limited ramjet applications. However, the advantages of 
boron slurry fuels are only available if combustion efficiencies comparable 
to those obtained with more conventional hydrocarbon fuels can be attained. 

In past combustor development work the realization of high levels of combustion 
efficiency has been achieved only at the cost of large total pressure losses 
in the combustor, which forces unacceptable design compromises. 

Because the mechanisms involved in boron slurry combustion are considerably 
different from those involved in the combustion of liquid hydrocarbon fuels, the 
empirical design guidelines that have been developed for conventional fuels 
do not apply. Even for conventional fuels, the need for reduced combustor 
volume and the packaging requirements that lead to the use of a sudden-expansion 
dump combustor have together reduced the applicability of empirical design 
and scaling rules. For this reason, an analytical tool for the design and 
scaling of sudden-expansion ramjet combustors, utilizing conventional hydro¬ 
carbon fuels has been developed and tested, as outlined in the preceding 
section. This code combines a semi-empirical analysis of the fuel injection 
processes with a sophisticated turbulent kinetic energy analysis of the mixing 
process and a quasiglobal finite-rate hydrocarbon chemical kinetics mechanism 
in an iterative analysis of a sudden-expansion combustor with large scale 
embedded recirculation zones. Because of current interest in boron slurry 
fuels for ramjets and the inapplicability of empirical guidelines for the 
design of ramjet combustors using these fuels, the development of a further 
modification to the modular model for liquid-fueled sudden-expansion ramjets, 
to extend the approach to boron fuel combustion, has been initiated under 
this AFOSR program. 

Boron slurry fuel combustion is a complex phenomenon involving two-phase 
flow and finite-rate reaction processes in the gas phase and on the particle 
surface. The boron is generally in the form of micron-size particles suspended 
in a carrier hydrocarbon fuel. One of the difficulties involved in boron 
combustion is that a surface oxide layer can form on the particle. This oxide 
has a boiling point below that of the metal, so that the boron combustion 
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process involves a diffusion of oxidizing gases through the oxide layer with 
a surface reaction mechanism involved at the particle surface. Since this is 
a relatively slow process, one strategy for obtaining high boron combustion 
efficiencies is to limit the amount of oxygen reaching the surface of the 
particle during the heat-up of the particle, i.e., carrying out the particle 
heat-up in a regime which is fuel-rich relative to the hydrocarbon carrier, 
and then exposing the hot particles and remaining fuel to additional air. 

This form of staged combustion is difficult to achieve in practice, and the 
assembly of a computer analysis to aid in achieving it is one of the goals 
of this task. 

The analysis of the combustion of a boron slurry fuel in a ramjet combustor 
begins just downstream of the point of injection of the fuel. A spray is 
formed, which involves a distribution of droplet sizes within which there is 
also a distribution of boron particle sizes. As the fuel evaporates and burns, 
there remains a cloud of boron particles, with varying sizes and varying oxide 
layer thicknesses. To compute the subsequent fate of these particles in the 
combustor, it can be assumed that the particles are in near dynamic equilibrium, 
which means that, with respect to the mean velocity in the combustor, the 
particles follow the gas phase, but they do not necessarily follow the turbulent 
fluctuations of the gas phase. Thus there is an effective diffusion of parti¬ 
cles relative to the gas phase. 

In the analysis, the continuous distribution of particles is discretized 
into a number of classes, with the class described by an average particle 
size and oxide layer thickness range. It is further assumed that these 
particles form a dilute continuum, so that the motion of each class can be 
described by a species transport equation and the thermal history of each 
class by a particulate energy equation. The coupling of these species and 
energy equations to the gas phase combustion and heat transfer mechanisms is 
through source terms which model the rate of consumption and heat transfer 
phenomena for each class of particles. Since the boron oxidation process is 
in general occurring simultaneously with the oxidation of the hydrocarbon 
carrier fuel, these chemical kinetic phenomena must be considered together, 
which requires the addition of species involved in the boron oxidation process 
to those already included in the model. 
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The two major modifications to the modular model that are required for 
the modeling of boron slurry combustion are the development and introduction 
of a particle-tracking methodology and the development and introduction of 
submodels for the boron particle heat-up and combustion processes. These two 
modifications are the subject of this Section. As there currently exists no 
comprehensive model for boron particle combustion, the two most complete 
models available have been coded as interchangeable subroutines which, in 
addition to their use within the modular model, can also be utilized as 
single-particle models. In the latter mode, these subroutines can be used 
in the development of more complete models of the boron combustion process, 
through comparison of single-particle model results with relevant experimental 
data.' The particle tracking formalism has been incorporated into the code 
utilizing a total of ten particle categories: as each additional category 
involves an additional species transport equation, in addition to an increase 
in the overall storage required for the code, ten categories have been selected 
as providing the best balance between the need for discrimination of particle 
size effects and the overall code run time required for solution. 

3.1 BORON IGNITION AND COMBUSTION MODELS 

There are a number of models for the ignition and combustion of single 
particles of boron. A variety of these have been investigated for possible 
application in the modular model for slurry-fueled ramjets, including those 
reported by King (Refs. 14, 15), Meese (Ref. 16), Vovchuck (Ref. 17), Mohan 
(Ref. 18) and Edelman (Refs. 19, 20). The component equations used in each 
of these models are summarized in Table 1, which also incorporates the 
assumptions built into the models; the nomenclature for Table 1 is given in 
Table 2. 

In reviewing these models, the emphasis has in each case been on the 
utility and generality of the approach with respect to the question of addressing 
the combustion of boron particles under the wide range of conditions appropriate 
to ramjet operation. The attention given to both diffusion controlled and 
kinetically controlled mechanisms was examined, as was the characterization 
(if any) of oxide removal by water. 
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In general, the formulations are based upon the conservation equations for 
mass and energy. The momentum equation is not explicitly coupled into the 
equation set but in certain of the models the relative velocity effect is 
accounted for through a Nusselt number correction on heat and mass transfer. 

In most of the formulations the energy equation is the central element of the 
model, with the size history and other transport relationships representing 
auxiliary elements of the model. 

King considers three intervals in the heat-up of the particle: below 
phase change of the boron, at phase change of the boron, and above the phase 
change of the boron (melting point 2450°K). Meese, who based his work on a 
1972 paper of King's, also includes the boil-off of the oxide layer. For 
conduction. King assumes the particle moves at the velocity of the surrounding 
gas, whereas Meese allows-for the incorporation in his calculation of the 
convective heat-transfer coefficient. Vovchuk includes no convective heat 
transfer and in Table 1, this section for Vovchuk is omitted. This is the 
form used when a section does not apply to an author’s work. Edelman also 
allows for convective heat transfer. The format for the listing of the 
radiation component in Table 1 is the same as that used for the convection terms. 
Each radiation equation could be interchanged with the other radiation equations 
with only minor modifications to match an author's model (i.e. volume based 
on a slab or a particle). It should be noted that Vovchuk includes a radiation 
term which he eventually neglects. 

For oxide formation, Kil® d V tfeese, andTfoharTuse the global reaction 
B + 3/402~-*~ 1 / 2 B 20 3 . Vovchuk allows for lower oxides BO, BO 2 , B 2 O 2 , forming 
at the surface. Edelman includes BO with B 2 0 3 but if BO forms it does not 
oxidize further. All reactions are diffusion controlled for the first four 
authors and kinetics controls the reactions Edelman uses. 

For oxide consumption, only King and Meese consider the chemical effect 
of water on boric oxide ( 82 ®x + H 20 -*~ 2 HB 02 ). The only other mechanism 
considered for oxide removal is evaporation. 

Diffusion equations in King's work are algebraic in form and depend upon 
absolute values of pressures at infinity and the surface. Vovchuk integrates 
the partial pressure from the surface to infinity. For diffusion through the 
oxide layer King and Meese have similar equations that vary only in one constant. 
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Meese's constant was derived from data unavailable to King and fits the data 
curve better. 

A major part of King's model is the phase change of the boron. Meese has 
a similar equation though derived in a slightly different fashion. Edelman 
assumes a finite evaporation rate of the boron. For evaporation of the oxide 
King and Meese are again similar. Vovchuk and Edelman have no oxide layer 
to remove. Mohan's phase change model, instead of being equilibrium controlled, 
is finite rate. 

Edelman is the only author who includes finite rate kinetics. He assumes 
the particle has ignited and studies surface oxidation kinetics. Vovchuk also 
studies only combustion in the diffusion controlled mode. 

The Nusselt number correlation allows for the effect of a velocity 
difference on convective heat transfer. King sets the value of two (N u = 2) by 
assuming a zero velocity difference. Meese and Edelman allow for variation in 
velocities. 

In Meese's model nothing else occurs while the oxide layer boils off 
leading to an independent expression for the oil off time. The end product of 
Vovchuk's work is the burning time. Mohan also derived an expression for the 
burning time. The radius and, in cases where the oxide is considered, the 
oxide layer thickness equation are basic to all but Vovchuk's model. 

Three of the models consider ignition: King, Meese, and Mohan. They all 
essentially take it to be the point where the oxide layer boils or flashes off. 
Other assumptions and restrictions made by each author appear at the end of 
Table 1. 

Based upon this review and formulation, the two most promising models 
for use in ramjet applications are those of King and Edelman. These are most 
general in their respective formulations; the others are either derivatives 
of them or are too restrictive to be usefully extended. 

3.2 BORON PARTICLE COMBUSTION SUBROUTINES FOR THE MODULAR MODEL 

The emphasis in this phase of the present work was on the creation of 
a single-particle combustion code which could employ either the methodology 
of King (Refs. 14, 15) or that employed by Edelman (Refs. 19, 20). The concept 
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TABLE 1. Summary of Available Boron Ignition and 
Combustion Models. 
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Heese 


The first stage of a two stage burn. During this time the oxide layer Is removed. 

King 

IGNITION 

A state at which no energy barriers stand between that state and steady-state combustion. 

I 

Only considers combustion. 

mm\ 


Similar to Xing's. 

Mohan 


N/A 

C^elran 


Oxide generation controlled by diffusion of oxygen through liquid oxide layer. 

Uniform particle temperature. 

Vaporitation and diffusion considered as series resistances. 

HpO removes B 2 O 3 by diffusion limited reaction (endothermic). 

Heat gain/loss by convection, radiation. 

Partial pressure of HBOp negligible in the free stream. 

Particle is stationary with respect to surrounding gasses. 

Paper only deals with ignition. 

King 

ASSUMPTIONS 

Initial B 2 O 3 lay®** thickness is uniform. 

Uniform particle temperature. 

The convective heat transfer coefficient, h, is known. 

Surroundings treated as black at temperature Traq. 

Prior to oxide removal, boron reaction rate is diffusion controlled and forms 62 O 3 . 

B 2 O 3 evaporation follows an Arrhenius type relation. 

Particle temperature remains constant during melting and boiling. 

Variation in h, surface reaction rate and oxide evaporation rate are neglected during oxide 
boiling. 

Water vapor reduces the activation energy of oxide evaporation. 

Initial particle temperature is below oxide boiling and boron melting. 

B 2 O 3 formed during boiling may be neglected. 

No net radiative exchange between particle and cloud. 

Convective heat transfer neglected In final computations. 

Heese 


The particle has been heated to equilibrium temperature. 

Reaction occurs at the surface. 

The diffusion coefficients for the oxidant and reaction products are Identical. 

Radiation can be neglected to a first approximation (d < 100yai). 

Particle temperature Is above oxide boiling and boron melting points. 

Vovchuk 


Diffusion controlled B ♦ 3/4 Op ■» 1/7 B^Oj. 
liquid layer Is Isothermal. 

Ideal gas. 

Absorptive equilibrium, for Op at C and ID steady-stage diffusion of Op within liquid layer. 

Plane two-dimensional slab model. 

Mohan 


No oxide layer. 

Uniform particle temperature. 

Surface combustion kinetics. 

Mnfte>rate evaporation of boron from the surface. 

Cdelsan 
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TABLE 2. Nomenclature for Expressions Summarized in Table 1 


2.1 King's Model 


TB1 

TBS 

tb 2 o 2 


liquid boron heat capacity (cal/gm°K) 
solid boron heat capacity (cal/gm°K) 
liquid boron oxide heat capacity (cal/gm°K) 


f 

h 

k 


(mw ) 8 

(mw)b 2 o 2 

Nu 


P 

Pr 


surf 

DO 

B 2°3 

P H 2 0, surf 
P HB0 2 , surf 
- 

q rx 

Q RX2 



fraction of boron in the liquid phase 

gas-particle heat transfer coefficient (cal/cm sec°K) 

mass transfer coefficient for transport of boric oxide gas from 

particle to free-stream (gm-mol/cm atm sec) 

boron atomic weight (gm/gm-mol) 

boric oxide molecular weight (gm/gm-mol) 

Nusselt Number 
total pressure (atm) 

oxygen partial pressure in free stream (atm) 

boric oxide partial pressure adjacent to particle surface (atm) 

boric oxide vapor pressure (atm) 

water gas partial pressure adjacent to particle surface (atm) 
HB0 2 partial pressure adjacent to particle surface (atm) 
boric oxide partial pressure in free-stream (atm) 
heat release of B(s) + |o 2 -*• ^B^O) (cal/gm-mol) 
heat release of B(l) + |o 2 -*• ^BgO^O) (cal/gm-mol) 

3 

gas law constant, 82.06 atm-cm /gm-mol°K 

molar rate of boron consumption (gm-mol/sec) 

molar evaporation rate of boric oxide (gm-mol/sec) 

molar rate of removal of by water reaction (gm-mol/sec) 

boron particle radius (cm) 

free stream gas temperature (°K) 

particle center temperature (°K) 

particle temperature (°K) 

surroundings radiation temperature (°K) 
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I 


^ 2^3 *^2 


u hbo 2 ,n 2 

d h 2 o,n 2 

AH., 


Pe 2°3 


particle surface temperature ('’K) 
molecular volume of species j (cm ) 
oxide layer thickness (cm) 

o 

diffusivity of gaseous boric oxide in nitrogen (cm /sec) 

2 

diffusivity of HB0 2 in nitrogen (cm /sec) 

2 

diffusivity of water gas in nitrogen (cm /sec) 

heat absorbed by reaction of H 2 0 with )(cal/gm-mol) 

heat of fusion of boron (cal/gm) 
heat of vaporization of B^^OMcal/gm-mol) 
evaporation coefficient of boric oxide liquid 
surroundings absorptivity 

2 

thermal diffusivity of boron (cm /sec) 

ip p 4 

Stefan-Boltzmann constant, 1.354*10 cal/cm sec°K 
particle emissivity 

Hertz-Knudsen impingement factor (gm-mol/cm atm sec) 
time (sec) 

3 

boron density (gm/crn ) 

3 

boric oxide density (gm/cm ) 


2.2 Model of Meese and Skifstad 

2 

A surface area, cm 


specific heat, cal/g-°K 

particle diameter, cm 

activation energy for evaporation 

heat of fusion of boron, cal/g 

heat of vaporization of boron oxide, cal/mole 

convective heat-transfer coefficient, cal/cm -sec-°K 

thermal conductivity, cal/cm-sec-°K 

molecular weight, g/g-mole 

mass, g 

pressure, atm 

Prandtl number of the gas 

heat of combustion of boron with diatomic oxygen to form liquid 


BgOj, cal/mole; total energy, cal 
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R 

r e 

Re 

r 

Ar 

T 

t 

x h 2 o 

a 

e 

X 

P 

a 


molar rate of consumption of boron, g-mole/sec; gas constant 
molar rate of evaporation of boron oxide, gm-mole/sec 
Reynolds number based on particle diameter and relative velocity 
between the gas and particle 
radius, cm 

oxide film thickness, r 2 -r.|, cm 
temperature, °K 
time, sec 

mole fraction of H 2 0 present in the gas 

surrounding absorptivity (=1.0) 
particle emissivity (=0.8) 
particle burning rate, cm/sec 
density, g/cm 

2 4 

Stefan-Boltzmann constant, cal/cm sec°K 


2.3 Mohan's Model 


C 

C t 

D 

d 

h 

k 

k* 

L 


average molecular velocity normal to the surface (cm/sec) 

heat capacity of B 2 0 3 (L)(cal/gm°K) 

diffusion coefficient (cm^/sec) 

diameter of particle (cm) 

thickness (cm) 

effective distribution coefficient related to the ratio of the 

mole fraction of 0 2 absorbed in B 2 0 3 (Jt) to gas phase mole fraction 

distribution ratio (liquid phase to gas phase 0, concentration) 

sum of heat losses by radiation and conduction — 2 - 

cfn sec 


P 

q B 

q C 

^B or C 
R° 

T 

t 

W 1 

X 


total ambient pressure (atm) 

heat liberated per unit mass of boron consumed at B (cal/g) 

heat absorbed per unit mass of B 2 0 3 vaporized at C (cal/g) 

regression rate of the surface (cm/sec) 

universal gas constant (8.316*10"^erg/gm-mole °K) 

temperature (°K) 

time (sec) 

atomic weight 

mole fraction of 0 2 in the ambient atmosphere 

ratio equilibrium vapor pressure of B 2 0 3 at Tc to total pressure 
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a evaporation coefficient 

e emissivity 

A gas phase thermal conductivity (cal/cnrsec °K) 

Pg q density boron and (gm/cm) 

t time for the new combustion regime to appear (sec) 

B,C surfaces 


2.4 


do 


AH. 


M„ 


P i 


R 

r 

T 

e 

A 


V i 


a 


T 


r 


Vovchuk's Model 

molar specific heat 

diffusion coefficient 

initial diameter of the boron droplet 

total enthalpy of species i 

heat of reaction of species i 

molar flux of species i 

mass flow of the boron reacting in unit time 
total pressure 

partial pressure of species i 

radiative heat flux 

universal gas constant 

radius 

temperature 

emissivity 

thermal conductivity of gas 
molecular weight 
stoichiometric coefficients 
density 

Stefan's constant 
droplet combustion time 


Subscripts 

B (molten) boron 

K particle or particle surface 

• infinity 


Note: units not stated, therefore any consistent set of units may be used. 






2.5 Edelman's Model 


A 

a 

B 

C r 


AH p 

K S 

L 

M„ 


Mr 


N2 


prod 


P 

P 

R 

T. 


V,B(T p ) 


w 


B,S 

f N 2 ,S 

^o 2 s 

r o 2 ,» 

Vrod.S 


a 

6 s 

e 

X 

Pr 


number of moles of oxygen per mole of product 
the total (dimensionless) particle consumption rate 
number of moles of boron per mole of product 
drag coefficient 

specific heat of the gas (ft^/sec^ °K) 

2 2 

specific heat of boron at the surface (ft /sec °K) 

heat of reaction (calorie/gram) 
surface reaction rate constant 
latent heat of vaporization 
molecular weight of boron (g/g-mole) 
molecular weight of nitrogen (g/g*mole) 
molecular weight of oxygen (g/g-mole) 

molecular weight of the product (g/g-mole) 
average molecular weight at the particle surface (g/g*mole) 
pressure (in atmospheres) 
vapor pressure of boron (atm) 

particle radius (feet) 
gas temperature (°K) 
particle temperature (°K) 
temperature of the wall (°K) 
gas velocity (ft/sec) 
particle velocity (ft/sec) 

mass fraction of boron at the surface of the particle 
mass fraction of nitrogen at the surface 

mass fraction of oxygen at the surface of the particle 

mass fraction of oxygen in the surrounding (gas) 

mass fraction of the product at the particle surface 
fraction of vaporized product 
particle bulk density (slugs/ft ) 
emissivity 

the thermal conductivity of nitrogen (lbs force/sec °K) 

o 

gas density (slugs/ft ) 
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3 

p density at the surface of the particle (slugs/ft ) 

j » 

a Stefan Boltzmann constant (lbs force/ft*sec °K j 
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was to develop a code which could stand alone as a tool to be used in the 
investigation of single-particle combustion phenomena, but which could also 
be easily incorporated into the modular model format without substantial 
modification. This dual capability was achieved by structuring the code as 
a series of subroutines which do not depend on the basic code in which they are 
found. In the stand-alone code they are called out of a main routine which 
does little more than obtain the needed inputs and call the first of them; in 
the modular approach this calling subprogram is a simple routine which corrects 
existing data to the units and variables required and calls the boron particle 
subroutine package. The subroutine package itself needs no modification in 
either case. 


Both models involve the solution of a set of coupled ordinary differential 
equations. The integration was carried out using subroutine HPCG, one of the 
standard differential equation solvers. HPCG uses Hamming's modified predictor- 
corrector method to obtain an approximate solution to the system. Although 
not self-starting (the initial solution is obtained using a Runge-Kutta 
technique) the technique is superior to the Runge-Kutta in that it requires 
evaluation of the right-hand side only twice per step. Furthermore, because 
the local truncation error falls out automatically calculation of the proper 
step size is facilitated (it has the capability of moving in either direction 
as the solution proceeds, but it never exceeds an input maximum). The 
technique used is a stable fourth order scheme. 

3.2.1 King's Model 

3.2.1.1 Model Description 

King's model includes the inhibiting effect of the boric oxide layer on 
ignition and explicitly calculates the particle melt time before ignition. 

It is thus the better model of the ignition process. However, because it 
does not include the kinetically controlled terms, it is poor for predicting 
results for ignited particles (indeed, it predicts a total thermal runaway). 

The basic model consists of four coupled equations relating particle 
temperature, particle size, oxide layer thickness and the fraction of the 
particle which is melted. It can be expressed as follows: 




*B Q RX ~ R E AH vap " R H aH H + 4 tt( V x)2 ^ h ^V T p^ +qca R^ T rad~ T p^ T <245Q f=( 


4 n r „ P n C + 4irr^X p R n C 
3 p P B p BS p B 9 0, p 


2 U 3 p B 2 0 3 


^ = 
dt 


0, T p = 2450, f<l 


(3) 


r b Qrx, - - R U AH U + 4ff(r n +x) 2 [h(VT n )+aca p (T; aH -Tj)] 


E vap ''H UM H 


i if p D C + 4nr 2 X p B n C 

p B P BL p B 2 0 3 p b ^ 0 ^ 


R rad P , T p >2450, f=l 


where 





Q RXZ 

r b 

r h 


T 


rad 


T 


X 


a R 

AH H 


liquid boric oxide heat capacity (cal/gm °K) 

liquid boron heat capacity (cal/gm °K) 
solid boron heat capacity (cal/gm °K) 

gas-particle heat transfer coefficient (cal/gm °K) 
heat release of B(S) + ^ ^ * 2 ^2®3^^ ca ^ gm ~ m0 ^ e ^ 

heat release of B(L) + f % 1 3 2 0 3 (£)(cal/gm-mole) 

molar rate of boron consumption (gm-mole/sec) 

molar rate of removal of B 2 0 3 by water reaction (gm-mole/sec) 

particle radius (cm) 

particle temperature (°K) 

surroundings radiation temperature (°K) 

free stream gas temperature (°K) 

oxide layer thickness (cm) 
surroundings absorptivity 

heat absorbed by the reaction of H 2 0 with B 2 0 3 U)(cal/gm-mole) 
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AH vap heat of vaporization of (cal/gm-mole) 

c particle emissivity 

O 

p b boron density (gm/cm ) 

3 

p R n boric oxide density (gm/cm ) 

b 2 u 3 

a Stefan-Boltzman Constant (cal/cm^sec °K 4 ) 

0, T p <2450, f=0 

R B Q RX~ R E AH vap'W 4 *^ V x ^ ^ h ( T °° ~ T p^ aea R [ T rad~ T p^2 . T p = 2450, f<l 

4/3^pP B AH m 

VO, T p >2450, f=l (4) 

where AH m = heat of fusion of boron (cal/gm), and the other variables are as 
defined above. 



^ 2 . = 
dt 


R b (mw) 


B 


4 " r p C B 


(5) 


where (MW) g = boron atomic weight (gm/gm-mole), and the other parameters are 
as above. Finally, 


dx 

dt 


(R b /2 - R e - R h )(MW) B203 

4,r p C B 2 0 3 


( 6 ) 


where (MW) D n is the boric oxide molecular weight (gm/gm-mole) and all other 
b 2 u 3 

variables are already defined. 

As can be seen, the King model is actually three models, of which one 
is selected based on current particle temperature, and melt history. This 
requires careful handling in the coding. The HPCG routine needs a subroutine 
to evaluate the derivatives at each time step. One cannot, however, merely 
wait until the particle temperature meets some criterion, then start using the 
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melt equation because the integrator may call this routine many times for every 
final solution it accepts. Furthermore, there is no easy way to tell while in 
the derivative routine whether or not the solution this time will be accepted. 
Thus it was necessary to write three essentially separate routines for the 
three states of the particle heating, melt, and burn. 

Initially the time range of interest is passed to the integrator, which 
starts with the assumption that the particle is not melted. When (and if) 
the particle reaches (or passes) the melt temperature (2450°K) the integration 
is stopped, values are set back to those at which Tp = 2450 exactly, and the 
integration procedure is initialized for the new problem (defined now as the 
calculation of the fraction of the boron which is liquid, f (Eq. 4)). When 
the liquid fraction reaches (or passes) unity, the integration is again stopped 
and the values are again reset (this time to those at which the liquid fraction 
is exactly 1.0) and the integration is initialized for the burn problem. 

This methodology allows automatic correct handling of those cases in which 
the particle never melts, those in which only partial melting occurs, as well 
as those cases in which ignition occurs. The implementation is all automatic 
and transparent to the user. 

Because the individual time steps are small with respect to the variables 
being calculated it was decided that a linear interpolation gave sufficient 
accuracy, and represented the best approach to the problem of resetting the 
variables at the boundaries of the particle states. A comparison of King's 
results with those obtained using the SAI code verifies the adequacy of this 
approach. 

The approach was checked by recalculating some of King's original results 
using the King model code. Two cases were selected for this on the basis of 
near-critical temperature determining whether or not ignition occurred. 

In both causes the particle radius is 10 microns and the oxide layer is 0.1 
microns. The initial particle temperature is 1800°K; the oxygen mole fraction 
is 0.2 and the pressure is 5 atmospheres. 

The first case has a gas temperature of 2100°K. Fig. 9 shows the 
particle temperature histories as calculated by King and by SAI. The results 
are essentially identical. The observed differences are more likely due to 
inaccuracies in reading the King curve from Ref. 16 than to be actual numerical 
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FIGURE 9. Temperature history for the 2100° case. 









differences. Fig. 10 compares the two calculations of the oxide layer thickness 
as a function of time. Again it is seen that the two calculations are essentially 
identical. Fig. 11 gives the particle radius history for the SAI calculation 
(King did not present this result). 

The second comparison case was run for an identical case except that the 
gas temperature was 2000°K. Here King shows that true ignition does not occur. 
Fig. 12 shows that the SAI code predicts essentially the same result as was 
obtained by King. Fig. 13 also confirms this result for the oxide layer 
thickness (the apparent difference again may well be due to difficulty in 
reading the detail on the plot of King's results - if the data were plotted 
in the same scale as they appear in Ref. 16 there would be no discernable 
difference). Comparison of Figs. 11 and 14 shows very little difference 
initially. However, in the degenerate combustion case the initial trend is 
continued whereas in the full ignition case (Fig. 11) the effect of particle 
heatup and melting is evident. 

The two cases taken together demonstrate that the King model has been 
faithfully reproduced in the SAI code. The system behaves just as predicted 
by King on both sides of the critical gas temperatures. 

3.2.2 Edelman's Model 

Edelman's model (Refs. 20, 21) neglects the inhibiting effect of the 
boric oxide layer on ignition and fails to calculate the particle melt time. 
However, it includes all the kinetically controlled mechanisms in the combustion 
equations. Consequently, it is superior to King's model for examination of 
what happens after the particle has ignited. 

The full Edelman model consists of nine coupled ordinary differential 
equations. However, as expanded by Edelman (Ref. 21) the system can be 
simplified to three equations without significant loss of accuracy; the 
remaining equations, although actually functions of time, can be treated 
algebraically. It is the simplified version which is here presented. The 
three differential equations are: 


ft. 3 wy v / 

dt ' 8R6 $ 


(7) 
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FIGURE 10. Oxide layer history for the 2100° case. 






KING MODEL RESULTS: R 






PARTICLE TEMPERATURE, 


2100 



FIGURE 12. Particle temperature history in the 2000° case. 
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FIGURE 13. Oxide layer thickness history for the 2000° case. 








FIGURE 14. Particle radius history in the 2000° case. 
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where = drag coefficient 

Pg = gas density (slugs per cubic foot) 
Vg = gas velocity (ft/sec) 

Vp = particle velocity (ft/sec) 

R = particle radius (feet) 

3 

and = particle bulk density (slugs/ft ) 


dR 

dt 


-axN 


2R6 C ’ 
s p 


( 8 ) 


where a = the total (dimensionless) particle consumption rate (see Eq. 10) 
x = the thermal conductivity of nitrogen (in pounds force/sec °K) 

Nu = Nusselt number 

2 2 

Cp = specific heat of the gas (in ft /sec °K), and the remaining 
variables as already defined. 


dt 


r 


3axNu 
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where 
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T 


P 


L 

A 




2 2 

specific heat of boron at the surface (ft /sec °K) 

gas temperature (°K) 

particle temperature (°K) 

latent heat of vaporization 

number of moles of oxygen per mole of product 

surface reaction rate constant 

o 

density at the surface of the particle (slugs/ft ) 
mass fraction of oxygen at the surface of the particle 

molecular weight of boron (g/g*mole) 

heat of reaction (calorie/gram) 
molecular weight of oxygen (g/g*mole) 
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o = Stefan-Boltzmann Constant (lb*force/ft*sec °K H ) 
c = emissivity 

T = temperature of the wall (°K) 

VV 

and other parameters are as defined above. 

The values for the remainder of the time-varying parameters are treated as 
simple algebraic equations (using the current value of Tp and R) as follows 


88.662C aR /mT M, 

_ 1 _ p _ v p 

XNu \ T v,B(T ) M b Y B,S K 


2M B K sV c C p P s RA 


M n XNu 
U 2 


where a = fraction of vaporized product 

Py j = vapor pressure of boron (atm) 

M g = average molecular weight at the particle surface (g/g-mole) 

Yg s = mass fraction of boron at the surface of the particle 
P = pressure (in atmospheres) 

and other variables are as previously defined. 


1 + BK. 


2,°° _ 

2(l-e~ a )C p R 


where Y n = mass fraction of oxygen in the surrounding (gas) 
>°° 

B = number of moles of boron per mole of product 
and the remaining variables are as previously defined. 
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where 


mass fraction of the product at the particle surface 
molecular weight of the product g/g*mole) 


Y 

prod.s 
^prod 

and the remaining parameters are defined above 


fl s = average molecular weight at the surface = 




+ 


Y 

prod,s J 

Vod 



where = molecular weight of nitrogen (g/g-mole) 

n 2 

Y n = mass fraction of nitrogen at the surface 


and the remaining variables as previously defined. 


(14) 



^ " Y 0 2 ,s ” Y B,s ~ Y prod,s^ 


(15) 


The Edelman model is, in a sense, more easily implemented than King’s 
because it does not treat separate realms of phenomena. However, as will be 
seen from an examination of Eqs. 10-13, virtually everything is a function of 
a, which itself is a function of both Y n . and Y 0 ... Hence one is faced with a 

Uo »S D > 5 

set of three highly non-linear equations^which are strongly coupled. 

In the original work done by Edelman an iterative solution technique was 
implemented which depended on the ratio of surface mass fractions to surround¬ 
ings mass fractions to converge to a solution. However, all attempts to 
implement this technique were met with uniform failure. Finally this approach 
was abandoned and a Newton-Raphson scheme was derived for the solution of 
the set of equations. This approach quickly led to convergence. However, a 
few troublesome anomalies appeared, the most serious of which was a calculated 

negative mass fraction for the product in some areas. Since this was extremely 
-9 

small (of order 10 ) it was decided to be due to inherent inaccuracy in the 

numerics (and/or machine roundoff) and zeroed out when it occurred in a final 
converged result (but not before). Consequently, the final convergence 
criteria were: a change in all three variables from one step to another of 
less than 10'^ or a change in all three variables of less than 0.1%. If either 






of these criteria were met, the solution was considered to have converged 

unless a negative mass fraction occurred. A significant negative mass fraction 

causes the routine to iterate again. (A significant negative is defined to 

-8 

be greater in absolute value than 10 ). 

Initial results from the Edelman model suggest that the diffusion terms 
are not the dominant ones in the equations. Hence the results suggest that 
the driver mechanism is kinetics; indeed, Fig. 15 shows that the particle 
radius decline is very close to linear with time. Even the temperature curve 
for the J7.3 micron particle is reasonably close to linear (Fig. 16). Results 
for a larger particle (50 micron diameter) are also nearly linear (Fig. 17), 
and a plot of the square of particle diameter vs. time (Fig. 18) does show some 
curvature. 

3.3 PARTICLE TRACKING FOR SLURRY FUEL COMBUSTION ANALYSES 

The combined effects of continuous reaction of boron particles and the 
non-uniform velocity and temperature fields within the combustor provide 
conditions resulting in a spectrum of particle sizes at every point in the 
flow. To completely specify the particle reaction rate and the production of 
gas phase species, the local particle size and nu.T.ber density distribution 
must be known. This in turn requires the development of a methodology for 
tracking the particle history throughout the flow. 

To track the particle history through the flow, the total amount of parti¬ 
culate material present at every point in the field is divided into a set of 
categories. These categories can be based on any parameter which can be used 
to differentiate particles: for the boron application both particle size 
and the thickness of the liquid layer (X) on the particle surface are 
relevant parameters. Both parameters may be used, with liquid layer thickness 
forming a subclassification in a basic particle size classification, but it 
should be kept in mind that each category used introduces a new variable in 
the overall combustor model with associated transport equation and storage 
requirements. 

For simplicity, consider a categorization based on particle size. Within 
each size range, all particles are combined to form an average with a volume 
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FIGURE IS. Particle radius history for the 17.3 p particle 
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FIGURE 17. Particle radius history for the 25 y particle. 








of particle diameter as a function of time. 











to surface area ratio equal to that of the summation over all contributions. 
This equivalent particle size is the Sauter mean value and is expressed as 


■ En 


,r f 3 /£n ir , 2 
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where n. particles have radius r.. The mass fraction of each category is 
equal to the sum of the mass fractions of each contribution. With the average 
size of category particles and the total category mass fraction known, the 
appropriate number of equivalent particles is established for every class. 

This representation provides an exact equivalence for the subsequent rate 
process provided the rate of change of mass depends directly upon the particle 
surface area. 

The particle distribution array in terms of size, mass fraction, and 
number in each category at every grid point can only be altered by changes 
in size of the particle due to chemical reaction (which can also, for boron, 
involve changes in liquid layer thickness) and by diffusion of particles from 
the neighboring flow. The radius of the particle changes in accordance with 
the particle consumption expression as the flowfield develops. Thus, after 
every step, the new particle size is given by 


r = r + 
o 



07) 


Where the step size. Ax, is evaluated based on the stability criteria incor¬ 
porated in the modular model. After every step the updated radius array 
must be reclassified in order to account for the transition of particles from 
one category to another. Since it is possible for more than one contribution 
to be made to a particular category, the total category mass must then be 
reaveraged using the Sauter criteria, to establish the proper radius and number 
for the group. The particle array in terms of mass fraction, size and number 
density resulting from the kinetics of the combustion process will be defined as 

G(J,I) k 

R(J,I) k at grid point 
N(J,I) k 




Each category in the particle array is altered by the diffusion of particles 
of the same category from neighboring upstream locations. Subdividing the 
particle mass fraction a. into each separate category and treating the category 
mass fractions as separate species, the diffusion equation can be written in 
finite difference form in von Mises coordinates as 


a i,n+l ,m 


,n,m 


+ 



Ax 

/( A *) 2 




(18) 


where the n,m subscripts refer to the grid point locations shown in Fig. 19. 


The second term in Eq. 18 represents the change in due to the kinetics 
process as the flow moves downstream from n to n+1. The sum of the first two 
terms is the total mass fraction of category particles present downstream due 
to kinetic processes, and the remaining term represents the contribution to 
the mass fraction from diffusion. Noting Fig. 19, the diffusional term can 
be split into the contributions from each of the three upstream points: 
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For each of the categories present at a given grid point, we may express the 
category mass fraction as 


6(0,I)n+l ,m 



+ aG£ + AG 3 


( 22 ) 
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where the kinetics contribution is evaluated as G(J,I)^ with its associated 
R(J,I)^ and N(J,I)| c . Knowing the category sizes associated with each upstream 
contribution the Sauter criteria (or other appropriate characterization 
criteria) can be used to define the new category radius: 

R(J,I)n+l,m = 


G(J,I)k + AG-j + AGg + AGg 
G(j ^Gj AG2 aG^ 

RptiTk + R(J,I)n,m+T + R(J,I)n,m + R(J,I)n,m-l 

Knowing the new mass fraction and size, the number density of particles is 
determined from: 


(23) 


N(J,I)n+l ,m 


pn+1,mG(J,I)n+l,m 
4/3up^R(J,I)n+l ,m 


(24) 


3.4 SUMMARY 

The particle tracking and boron particle combustion subroutines described 
in this section represent major developments with respect to the assembly of a 
modular model formulation for the boron-slurry fueled ramjet combustor. 

However, considerable further work still remains to complete the modular model 
development. Some of this work involves the definition of improved models 
for the kinetics of the boron combustion process and thus requires additional 
experimental data. 

As was noted in the discussion of boron particle combustion models in 
this section, neither of the models presently available is complete: King's 
model (Refs. 14, 15) ignores chemical kinetic effects and describes only the 
ignition stage, while Edelman's model (Refs. 19, 20), which involves kinetically 
controlled combustion phenomena, ignores the oxide layer in its formulation. 

Of the two, the most easily extended is the Edelman model, and the generaliza¬ 
tion of this model to incorporate the effects of the oxide layer on the combus¬ 
tion process is presently being investigated. 

Further development with respect to the modeling of two-phase flow 
phenomena within the context of the modular formulation involves the 
incorporation within the model of the effects of thermal nonequilibrium and 
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modeling of the simultaneous oxidation of boron-containing species and the 
hydrocarbon species in the gas phase. Work in these areas is also in progress. 








4. FORMULATION OF A MATHEMATICAL MODEL OF A DUCTED ROCKET COMBUSTOR 


The ducted rocket has certain features that are closely related to 
ramjet systems, particularly integral rocket ramjets which utilize the 
sudden-expansion ("dump") combustor concept coupled with multiple side entry 
air inlets. The ducted rocket is sometimes referred to as a gas generator 
or solid propellant rocket-fueled ramjet. Actually, it is the rocket feature 
that constitutes the major differentiating feature between ducted rockets 
and liquid fueled ramjets, since the rocket exhaust, which provides the fuel 
injection into the ramburner, is choked and its momentum is significant. 

Because of the underexpansion of the rocket exhaust local supersonic flow 
is encountered and the significant level of the rocket exhaust momentum 
relative to the airflow momentum can strongly affect the overall ramburner 
flowfield. These factors are not present in a liquid fueled ramjet. Signifi¬ 
cant differences also exist between the ducted rocket fuel, which is made up 
of the products of combustion of a fuel-rich solid propellant, and the liquid 
(or slurry) fuels used in ramjets. Further, while practical sudden expansion 
ramjet combustors can be axisymmetric, ducted rocket configurations, as a 
rule, are three-dimensional. Thus the analytical simplifications offered 
by axisymmetry are not available in modeling the ducted rocket flowfield. 

Therefore, although there is a logical connection between the combustion 
processes in ramjets and in ducted rockets, it is clear that there are 
essential differences that must be addressed in order to establish a better 
quantitative understanding of the mechanisms and parameters involved in the 
ducted rocket combustion process. Despite the evident difficulties, develop¬ 
ment problems that have been encountered clearly show the need for research 
into ducted rocket combustor aerodynamics: flow patterns, mixing rates, 
and chemical kinetic rates in ducted rocket environments need to be defined, 
and analyses tailored to the specific characteristics of the ducted rocket 
combustion chamber need to be developed. 

The complexity inherent in the ducted rocket combustion chamber argues 
for the use of a modular formulation such as has been successfully developed 
for axisymmetric liquid fueled sudden expansion ramjet combustor configurations. 
However, while a data base which provides some description of the combustor 


59 



aerodynamics exists for the dump combustor, no such data base exists in the 
case of the ducted rocket. Since the modular approach depends on an a priori 
description of the size and location of major flowfield elements, particularly 
regions of recirculating flow, an initial data base is critical to the 
development of a detailed modular model. In the absence of experimental 
data, some information on the aerodynamics involved in the ducted rocket 
combustion chamber flowfield can be obtained through the use of a three- 
dimensional Navier-Stokes solution procedure. However, this type of modeling 
is complex and not readily amenable to the incorporation of the detailed 
chemical kinetics which is also required to gain an overall analytical descrip¬ 
tion of ducted rocket flowfield details. It is, on the other hand, possible 
to use the Navier-Stokes solution procedure either for nonreacting flow or, 
using a simple representation of the primary heat release reactions, for 
reacting flow, to provide the detailed information necessary to construct an 
appropriate modular model. This is the approach adopted in the present work. 

In order to discuss the requirements on a three-dimensional Navier-Stokes 
solution procedure for the ducted rocket it is useful to consider a specific 
configuration. To this end, the configuration which is to be utilized in an 
experimental program at the Aeropropulsion Laboratory (AFWAL) is of interest: 
a sketch of this configuration is shown in Fig. ZO. Basic modeling requirements 
for the main combustor include the use of the three-dimensional form of the 
Navier-Stokes equations in cylindrical coordinates, along with a suitable 
turbulence model and the capability to provide the particle tracking required 
to compute two-phase flow transport phenomena as outlined in the previous 
section. For this specific configuration a half-domain computation is 
sufficient but for generality the capability to compute the full flowfield 
domain may be necessary. 

Two significant complications are inherent with this flowfield. First, 
the air inlet to the combustor is provided through rectangular ducts, so that 
the inlet airflow is quite nonuniform, and second, the flowfield in the region 
of the gas generator nozzles must also be computed. In the latter case, the 
scale of the flow is significantly smaller than the scale which characterizes 
the overall flowfield. To adequately handle both of these complications, 
variable spacing capability is required using both nonuniform mesh spacing and 
(possibly) coordinate transformations. 
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FIGURE 20. Ducted rocket combustor configuration, 
AFWAL/APL. 
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The flowfield within the inlet ducts involves, since recirculations are 
not expected to be present, a three-dimensional parabolic or quasiparabolic 
flow. In general, the flow in this duct would be coupled to the conditions 
within the combustor, but at least as a first approximation an uncoupled 
analysis can be used. One aspect of this flowfield that requires special 
attention is the matching of the inlet duct grid with that utilized for the 
main combustor, as shown schematically in Fig. 21. Matching of grids is 
also a requirement with respect to the injector nozzle face, as also shown on 
Fig. 21. The most suitable approach to this aspect of the problem, given the 
different scales which characterize the injector nozzle and the combustor as 
a whole is the use of a control volume analysis for the computational region 
surrounding the injector nozzles. 

Development of an aerodynamic model for the combustor shown schematically 
in Fig. 20 has been based, initially, on the use of the three-dimensional 
elliptic formulation developed at Imperial College, and described, for example, 
by Patankar (Ref. 21). This code is a primitive-variables approach, which is 
to say that it solves the equations for the mean velocity components directly, 
rather than using transformed variables such as the stream function and 
vorticity. The solution is obtained iteratively using an underrelaxation 
technique, with a semi-implicit pressure-correction procedure and a velocity 
update scheme. For the ducted rocket application the code requires considerable 
modification, and during this modification effort the structure of the 
technique is being revised to allow modifications to various aspects of the 
model (i.e., turbulence models, two-phase flow formulations, and simple 
chemical kinetics approaches) to be more readily incorporated. The basic code 
is, however, a steady state formulation; an unsteady-state approach may 
eventually be required because of the oscillatory nature of the flow that has 
been observed in some ducted rocket combustor experiments. 
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FIGURE 21. Computational grid requirement for ducted rocket 
combustor computation. 









5. OVERALL SUMMARY AND FUTURE WORK 


Research effort for this period has focused on three major areas of 
investigation: further examination of the use of the modular ramjet combustor 
model as a technique for the analysis and interpretation of combustor test 
data, the development and implementation of a modular model for the boron 
slurry fueled ramjet, and the development of analytical models for ducted 
rocket combustors. Further examination of the liquid-fueled ramjet modular 
model has shown that the stirred-reactor representation of the flame-stabilizing 
recirculation zones in a dump combustor provides a useful and powerful tool 
for the examination of dump combustor flame stability characteristics. 

This work has also demonstrated the utility of the model in terms of the 
analysis and interpretation of ramjet test data. One area of required future 
work has been shown to involve the engineering modeling of chemical kinetics 
phenomena under fuel-rich conditions: in some ramjet designs and under some 
operating conditions local regions within the combustor, and particularly 
the recirculation regions, can be considerably more fuel-rich than would be 
expected based on the overall fuel-air ratio for the operating condition. 
Engineering chemical kinetics models for fuel-rich conditions are being 
developed by SAI under other research programs, and future work in this area 
of the current AFOSR program will involve the incorporation of these models 
into the modular formulation of the liquid-fueled sudden-expansion ramjet. 

Initial work relevant to the boron-slurry fueled ramjet has centered 
on the analysis and implementation of existing models for boron particle 
combustion and the development and implementation of a particle tracking 
technique for use with the modular model formulation. Existing models for 
boron particle ignition and combustion do not encompass the full range of 
combustion phenomena relevant to boron combustion necessary to provide a 
comprehensive model. Future work in this area involves the extension of 
existing approaches to consider both diffusion- and chemical kinetically- 
controlled combustion phenomena, including the effects of an oxide surface 
layer on the particle. The particle tracking methodology developed for use in 
the modular model provides the only realistic method of computing the combus¬ 
tion behavior of particle clouds, but this approach needs to be extended to 
include thermal nonequilibrium between the particles and gas phase and the 
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simultaneous combustion of boron and gas-phase hydrocarbon species. Future 
work will involve these extensions. 

An aerodynamic model for the ducted rocket combustor configuration is 
under development, based initially on the three-dimensional elliptic formula¬ 
tion developed at Imperial College. Computations of the flowfield within a 
representative ducted rocket configuration will be carried out during the 
next year's work. The results of these computations will be used to assess 
the utility of the approach (including the possibility of there being a 
requirement for the development of a time-dependent model) and to provide a 
baseline of flowfield detail required to initiate the development of a detailed 
modular model of this flowfield. 


65 


F 









6. REFERENCES 


1. Harsha, P. T., and Edelman, R. B., "Mixing and Combustion in High Speed 
Air Flows," AFOSR TR 80-1005, SAI-80-020-CP, Science Applications, Inc., 
Canoga Park, CA, July 1980. 

2. Edelman, R. B., and Fortune, 0., "A Quasiglobal Chemical Kinetic Model 
for the Finite-Rate Combustion of Hydrocarbon Fuels," AIAA Paper 69-86, 

1969. 

3. Launder, B. E., Morse, A., Rodi, W., and Spalding, D. B., "Prediction of 
Free Shear Flows - A Comparison of the Performance of Six Turbulence 
Models," Free Turbulent Shear Flows, Vol. I, Conference Proceedings , NASA 
SP-321 , 1973, pp." 463-519^ 

4. Edelman, R. B., Harsha, P. T., and Schmotolocha, S. N., "Modeling 
Techniques for the Analysis of Ramjet Combustion Processes," AIAA Paper 
80-1190, AIAA/SAE/ASME 16th Joint Propulsion Conference, Hartford, 

Conn., June 30-July 2, 1980. Also AIAA Journal , Vol. 19, No. 5, May 1981, 
pp. 601-609. 

5. Harsha, P. T., and Edelman, R. B., "AFOSR Interim Scientific Report: 

Mixing and Combustion in High Speed Air Flows," A0FSR-TR-79-0912, 
SAI-79-013WH, Science Applications, Inc., Canoga Park, CA, April 1979. 

6. Harsha, P. T., and Edelman, R. B., "Application of Modular Modeling to 
Ramjet Performance Prediction," AIAA Paper 78-944, 14th Joint Propul¬ 
sion Conference, Las Vegas, Nev., July 25-27, 1978. 

7. Curran, E. T., "An Investigation of Flame Stability in a Coaxial Dump 
Combustor," PhD. Thesis, Air Force Institute of Technology, 1979. 

8. Schmotolocha, S. N., and Economos, C., "An Experimental Combustion and 
Flame Stabilization Study of Dump Burners," AFOSR-TR-75-1446, Sept. 1975. 

9. Rosfjord, T. J., "Experimental Investigation of Fuel Sprays Formed by 
Ramburner Injectors," 1979 JANNAF Propulsion Meeting, Vol. 2, CPIA Pub. 300, 
1979, pp. 585-607. 

10. Pelmas, R. J., et al., "Integral Boost Ramburner Technology for Volume- 
Limited Missiles," Vol. I, AFAPL-TR-73-55, July 1973. 

11. Edelman, R. B., and Harsha, P. T., "Laminar and Turbulent Gas Dynamics 

in Combustors - Current Status," Progress in Energy and Combustion Science, 
Vol. 4, No. 1, 1978, pp. 1-62. 

12. Hautman, D. J., Dryer, F. L., Schug, K. P., and Glassman, I., "A Multistep 
Overall Kinetic Mechanism for the Oxidation of Hydrocarbons," Presented 

at 1980 Technical Meeting, Eastern Section, The Combustion Institute, 

Nov. 12-14, 1980. 




13. "Fundamental Characterization of Alternative Fuel Effects in Continuous 
Combustion Systems," SAI QPR 029-1, April 15, 1981. 

14. King, M. K., "Boron Particle in Hot Gas Streams," Combustion Science and 
Technology , Vol. 8, 1974, pp. 255-273. 

15. King, M. K., "Boron Ignition and Combustion in Air-Augmented Rocket After¬ 
burners," Combustion Science and Technology , Vol. 5, 1972, pp. 155-164. 

16. Meese, R. A. and Skifstad, J. G., "Ignition and Global Combustion Models 
for Clouds of Boron Particles," AIAA Journal , Vol. 12, No. 1 (1974), 

pp. 71-77. 

17. Vovchuk, Y. I. et al., "High Temperature Combustion of an Immobile Boron 
Particle in an Oxygen-Bearing Medium," Combustion, Explosion and Shock 
Waves , Vol. 11 (1975), pp. 471-476. 

18. Mohan, G. and Williams, F. A., "Ignition and Combustion of Boron in 0 J 
Inert Atmospheres," AIAA Journal , Vol. 10, No. 6 (1972), pp. 776-783/ 

19. Edelman, R. B. et al., "Part III: Theoretical and Experimental Investiga¬ 
tion of Metal Particle Combustion in a Supersonic Gas Stream," Gas Genera¬ 
tor Fueled Scramjet Program , 1971, AFAPL-TR-70-84, pp. 20-37. 

20. Edelman, R. B., Genovese, J., and Fortune, 0. F., "Some Aspects of Two- 
Phase Flows with Mixing and Combustion in Bounded and Unbounded Flows," 
Journal of Spacecraft and Rockets , Vol. 8, No. 4 (1971), pp. 82-87. 

21. Patankar, S. V., "Numerical Prediction of Three-Dimensional Flows," 

Studies in Convection , Vol. 1, B. E. Launder, ed.. Academic Press, New York, 
1975, pp. 1-79. 


67 






